function stats = CalcCrossSecMoments(parms, prob_Nxz, Phi_Nxz, stats_in)

    stats = stats_in;
    
    Phi_Nxz_1y = mpower2(Phi_Nxz,12);
    Phi_Nxz_5y = mpower2(Phi_Nxz_1y,5);

    x_Nxz = permute(repmat(parms.grid_x,[1 parms.NN parms.Nz]),[2 1 3]);
    U_Nxz = repmat(1 - parms.grid_N(:),[1 parms.Nx parms.Nz]);
    REC_Nxz = NaN*zeros(parms.NN,parms.Nx,parms.Nz);
    if parms.Nz==2
        REC_Nxz(:,:,1) = 1;
        REC_Nxz(:,:,2) = 0;
    else
        % N state model
        REC_Nxz(:,:,1:parms.IDX_RECESSION_THRESH) = 1;
        REC_Nxz(:,:,parms.IDX_RECESSION_THRESH+1:end) = 0;
    end
    
    phi_of_x = @(x) 1./(1 + exp(-x));
    phi_Nxz = phi_of_x(x_Nxz);
    log_phi_Nxz = log(phi_Nxz);
    
    [mean_phi,std_phi] = Calc_monthly_vol(prob_Nxz,phi_Nxz);
    
    stats.cross_sec_cons.mean_phi = mean_phi;
    stats.cross_sec_cons.std_phi = std_phi;
    
    mean_phi_given_z = zeros(1,parms.Nz);
    
    for iz = 1:parms.Nz
        prob_x_given_z = sum(prob_Nxz(:,:,iz),1);
        prob_x_given_z(:) = prob_x_given_z(:)/sum(prob_x_given_z);
        mean_phi_given_z(iz) = sum(prob_x_given_z(:).*phi_of_x(parms.grid_x(:)));
    end
    
    stats.cross_sec_cons.mean_phi_given_z = mean_phi_given_z;
    
    % percentiles of phi
    cdf_x = cumsum(squeeze(sum(prob_Nxz,[1 3])));
    cdf_x = cdf_x(:)/cdf_x(end);
    
    Pvals = [5 25 50 75 95]/100;
    pctiles = zeros(size(Pvals));
    
    for i = 1:numel(Pvals)
        idx = find(cdf_x(1:end-1)<=Pvals(i) & Pvals(i)<=cdf_x(2:end),1,'first');
        if cdf_x(idx)>=cdf_x(idx+1)
            % just pick the mid point
            pctiles(i) = phi_of_x((parms.grid_x(idx) + parms.grid_x(idx+1))/2);
        else
            pctiles(i) = phi_of_x(interp1(cdf_x([idx; idx+1]), parms.grid_x([idx; idx+1]), ...
                Pvals(i)));
        end
    end
    
    stats.cross_sec_cons.phi_Pvals = Pvals;
    stats.cross_sec_cons.phi_pctiles = pctiles;
    
    temp_Nxz_for_mu2 = (1 - U_Nxz).*U_Nxz.*(log_phi_Nxz.^2);
    [stats.cross_sec_cons.mean_mu2_1y,stats.cross_sec_cons.std_mu2_1y] = ...
        CalcTimeSeriesSumVol(temp_Nxz_for_mu2,prob_Nxz,Phi_Nxz,12,Phi_Nxz_1y);
    [stats.cross_sec_cons.mean_mu2_5y,stats.cross_sec_cons.std_mu2_5y] = ...
        CalcTimeSeriesSumVol(temp_Nxz_for_mu2,prob_Nxz,Phi_Nxz,60,Phi_Nxz_5y);

    temp_Nxz_for_mu3 = (1 - U_Nxz).*U_Nxz.*(1 - 2*U_Nxz).*(log_phi_Nxz.^3);
    stats.cross_sec_cons.std_mu3_1y = CalcTimeSeriesDiffVol(temp_Nxz_for_mu3,prob_Nxz,Phi_Nxz,12,Phi_Nxz_1y);
    stats.cross_sec_cons.std_mu3_5y = CalcTimeSeriesDiffVol(temp_Nxz_for_mu3,prob_Nxz,Phi_Nxz,60,Phi_Nxz_5y);
    
    % autocorrelation and autocovariance of mu3
    % 1. corr(mu3(t->t+12),mu3(t+12->t+12n)), n = 2,3,...,6
    % 2. cov(mu3(t->t+12),mu3(t+12->t+12n)), n = 2,3,...,6
    % 3. conditional version of 1 and 2
    
    corr_mu3_1y_ny = zeros(1,5);
    cov_mu3_1y_ny = zeros(1,5);
    
    if parms.Nz==2
        % 2 state model
        corr_mu3_1y_ny_given_z = zeros(parms.Nz,5);
        cov_mu3_1y_ny_given_z = zeros(parms.Nz,5);
    else
        % N state model
        corr_mu3_1y_ny_bus_cycle = zeros(2,5);
        cov_mu3_1y_ny_bus_cycle = zeros(2,5);
    end
    
    % correlation between mu3 and change in U
    % 1. corr(mu3(t->t+12n),U(t+12n)-U(t)) for n=1,2,...,5
    % 2. cov(mu3(t->t+12n),U(t+12n)-U(t)) for n=1,2,...,5
    
    corr_mu3_dU_ny = zeros(1,5);
    cov_mu3_dU_ny = zeros(1,5);
    
    % correlation between mu3 and recessions
    % 1. corr(mu3(t->t+12n),REC(t+12n)) for n=1,2,...,5
    % 2. cov(mu3(t->t+12n),REC(t+12n)) for n=1,2,...,5
    
    corr_mu3_REC_ny = zeros(1,5);
    cov_mu3_REC_ny = zeros(1,5);
    
    Phi_Nxz_n = Phi_Nxz_1y;
    for n = 1:5
        [corr_mu3_1y_ny(n),cov_mu3_1y_ny(n),~,~] = CalcTimeSeriesCorr1(temp_Nxz_for_mu3,prob_Nxz,Phi_Nxz,12,12*n,...
            Phi_Nxz_1y,Phi_Nxz_n);
        if parms.Nz==2
            % 2 state model
            [corr_mu3_1y_ny_given_z(:,n),cov_mu3_1y_ny_given_z(:,n)] = CalcTimeSeriesCorr1_cond_on_z(temp_Nxz_for_mu3,prob_Nxz,Phi_Nxz,12,12*n,...
                Phi_Nxz_1y,Phi_Nxz_n);
        else
            % N state model
            [cov_mu3_1y_ny_bus_cycle(:,n),corr_mu3_1y_ny_bus_cycle(:,n)] = ...
                CalcTimeSeriesCorr1_cond_on_bus_cycle(parms.IDX_RECESSION_THRESH,temp_Nxz_for_mu3,prob_Nxz,Phi_Nxz,12,12*n,...
                Phi_Nxz_1y,Phi_Nxz_n);
        end
        [corr_mu3_dU_ny(n),cov_mu3_dU_ny(n),~,~] = CalcTimeSeriesCorr2(temp_Nxz_for_mu3,U_Nxz,prob_Nxz,Phi_Nxz,12*n,Phi_Nxz_n);
        [corr_mu3_REC_ny(n),cov_mu3_REC_ny(n),~,~] = CalcTimeSeriesCorr3(temp_Nxz_for_mu3,REC_Nxz,prob_Nxz,Phi_Nxz,12*n,Phi_Nxz_n);
        
        Phi_Nxz_n = Phi_Nxz_n*Phi_Nxz_1y; % update
    end
    
    stats.cross_sec_cons.corr_mu3_1y_ny = corr_mu3_1y_ny;
    stats.cross_sec_cons.cov_mu3_1y_ny = cov_mu3_1y_ny;
    
    if parms.Nz==2
        % 2 state model
        for iz = 1:parms.Nz
            stats.cross_sec_cons.(['corr_mu3_1y_ny_given_z',num2str(iz)]) = corr_mu3_1y_ny_given_z(iz,:);
            stats.cross_sec_cons.(['cov_mu3_1y_ny_given_z',num2str(iz)]) = cov_mu3_1y_ny_given_z(iz,:);
        end
    else
        % N state model
        stats.cross_sec_cons.corr_mu3_1y_ny_given_recession = corr_mu3_1y_ny_bus_cycle(1,:);
        stats.cross_sec_cons.corr_mu3_1y_ny_given_expansion = corr_mu3_1y_ny_bus_cycle(2,:);
        stats.cross_sec_cons.cov_mu3_1y_ny_given_recession = cov_mu3_1y_ny_bus_cycle(1,:);
        stats.cross_sec_cons.cov_mu3_1y_ny_given_expansion = cov_mu3_1y_ny_bus_cycle(2,:);
    end
    
    stats.cross_sec_cons.corr_mu3_dU_ny = corr_mu3_dU_ny;
    stats.cross_sec_cons.cov_mu3_dU_ny = cov_mu3_dU_ny;
    
    stats.cross_sec_cons.corr_mu3_REC_ny = corr_mu3_REC_ny;
    stats.cross_sec_cons.cov_mu3_REC_ny = cov_mu3_REC_ny;

end